Inverse Optimization Techniques for Targeted Self- Assembly 
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This article reviews recent inverse statistical-mechanical methodologies that we have devised to 
optimize interaction potentials in soft matter systems that correspond to stable "target" structures. 
We are interested in finding the interaction potential, not necessarily pairwise additive or spherically 
symmetric, that stabilizes a targeted many-body system by generally incorporating complete con- 
figurational information. Unlike previous work, our primary interest is in the possible many-body 
structures that may be generated, some of which may include interesting but known structures, 
while others may represent entirely new structural motifs. Soft matter systems, such as colloids 
and polymers, offer a versatile means of realizing the optimized interactions. It is shown that these 
inverse approaches hold great promise for controlling self-assembly to a degree that surpasses the 
less-than-optimal path that nature has provided. Indeed, we envision being able to "tailor" poten- 
tials that produce varying degrees of disorder, thus extending the traditional idea of self-assembly 
to incorporate both amorphous and crystalline structures as well as quasicrystals. The notion of 
tailoring potentials that correspond to targeted structures is motivated by the rich fundamental 
statistical-mechanical issues and questions offered by this fascinating inverse problem as well as our 
recent ability to identify structures that have optimal bulk properties or desirable performance char- 
acteristics. Recent results have already led to a deeper basic understanding of the mathematical 
relationship between the collective structural behavior of many-body systems and their interac- 
tions, as well as optimized potentials that enable self-assembly of ordered and disordered particle 
configurations with novel structural and bulk properties. 



I. INTRODUCTION 

The term "self-assembly" typically describes pro- 
cesses in which entities (atoms, molecules, aggregates of 
molecules, etc.) spontaneously arrange themselves into a 
larger ordered and functioning structure. Biology offers 
wonderful examples, including the spontaneous forma- 
tion of the DNA double helix from two complementary 
oligonucleotide chains, the formation of lipid bilayers to 
produce membranes, and the folding of proteins into a 
biologically active state. 

On the synthetic side, molecular self-assembly is a po- 
tentially powerful method to fabricate nanostructures as 
an alternative to nanolithography. For example, it has 
been demonstrated that intricate two-dimensional struc- 
tures can emerge by the placement of organic molecules 
onto inorganic surfaces. 1 Block copolymers can self- 
assemble into ordered arrays that have possible use as 
photonic band-gap materials 3 Self-assembly based on 
contact electrification seems to be a powerful means to or- 
ganize macroscopic dielectric particles of various shapes 
into extended, ordered structures^ Highly robust self- 
assembly of unique, small clusters of microspheres that 
can themselves be used for self-assembly of more com- 
plex architectures has been demonstrated^ It has been 
shown that gold nanowires can be assembled by func- 
tionalizing nanoparticles with organic molecules^ DNA- 
mediated assembly of micrometer-size polystyrene parti- 
cles in solution could enable one to build complex struc- 
tures starting from a mesoscale template or seed fol- 



lowed by self-assembly^ These examples offer glimpses 
into the materials science of the future - devising build- 
ing blocks with specific interactions that can self-organize 
on a larger set of length scales. 

This is an emerging field with a wealth of experimen- 
tal data that has been supported theoretically and com- 
putationally using the "forward" approach of statisti- 
cal mechanics ^^ i 10 ' 1:L i 12 ' 13 i 14 i 15 ' 16 Such an approach has 
generated a long and insightful tradition. The forward 
approach identifies a known material system that pos- 
sesses scientific and/or technological interest, creates a 
manageable approximation to the interparticle interac- 
tions that operate in that material, and exploits simula- 
tion and analytical methods to predict non-obvious de- 
tails concerning structural, thermodynamic and kinetic 
features of the system. 

In the last several years, inverse statistical-mechanical 
methods have been devised that find optimized interac- 
tions that most robustly and spontaneously lead to a 
targeted many-particle configuration of the system for a 
wide range of conditions, 17 ' 18 ! 19 ! 20 ! 21 ! 22 ! 23 ! 24 ! 25 ! 26 ' 27 This 
article reviews these nascent developments as well as 
other closely related inverse realizability problems that 
we have introduced^^i^ 3 ^ 4 .^^ 7 . all of which 
are solved using various optimization techniques. 

Results produced by these inverse approaches have al- 
ready led to a deeper fundamental understanding of the 
mathematical relationship between the collective struc- 
tural behavior of many-body systems and the interac- 
tions: a basic problem in materials science and con- 
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densed matter theory. As will be shown, such methodolo- 
gies hold great promise to control self-assembly in many- 
particle systems to a degree that surpasses the less-than- 
optimal path that nature has provided. Indeed, employ- 
ing such inverse optimization methods, we envision being 
able to "tailor" potentials that produce varying degrees 
of disorder, thus extending the traditional idea of self- 
assembly to incorporate both amorphous and crystalline 
structures as well as quasicrystals. 

Output from these optimization techniques could then 
be applied to create de novo colloidal particles or polymer 
systems with interactions that yield these structures at 
the nanoscopic and microscopic length scales. Thus, this 
work has important implications for the future synthesis 
of novel materials. 

Colloidal particles suspended in solution provide an 
ideal experimental testbed to realize the optimized po- 
tentials, since both repulsive and attractive interac- 
tions can be tuned (e.g., via particle surface modifi- 
cation or the addition of electrolytes) ^±3SJ9joj^42 
and therefore offer a panoply of possible potentials 
that far extends the range offered by molecular sys- 
tems. Effective pair interactions in colloids can contain 
hard-core, charge dispersion (van der Waals), dipole- 
dipole (electric- or magnetic- field induced^ 2 -), screened- 
Coulombic (Yukawa), and short-ranged attractive deple- 
tion contributions. Polymer systems also offer a versatile 
means of realizing optimized soft interactions ! 11 ^ 3 

The idea of tailoring potentials to generate targeted 
structures is motivated by the rich array of fundamen- 
tal issues and questions offered by this fascinating in- 
verse statistical-mechanical problem as well as our recent 
ability to identify the structures that have optimal bulk 
properties or desirable performance characteristics. The 
latter includes novel crysta l 44 ' 45 ' 46 and quasicrysta l 25 ' 47 
structures for photonic band-gap applications, ma- 
terials with negative or vanishing thermal expan- 
sion coefficients ; 48 ' 49 ' 50 materials with negative Poisson 
ratios ; 26 ' 51 ' 52 ' 53 ' 54 ' 55 ' 56 materials with optimal or novel 
transport and mechanical properties ; 57 ' 58 ' 59 ' 60 ' 61 ' 62 ' 63 
mesoporous solids for applications in catalysis, separa- 
tions, sensors and electronics ^J 4 ^ and systems character- 
ized by entropically driven inverse freezing ; 66 ' 67 to men- 
tion a few examples. 

The recent inverse techniques that are reviewed 
here differ from so-called "reverse" Monte Carlo 
method a 68 ' 69 ' 70 ' 71 in several important respects. The lat- 
ter techniques are concerned almost invariably with ob- 
taining a spherically symmetric pair potential from an 
experimentally observed structure factor (as measured 
from scattering experiments) or real-space pair correla- 
tion function usually for stable liquid phases or glassy 
states of matter. By contrast, our interest is in finding 
the interaction potential, not necessarily pairwise addi- 
tive or spherically symmetric, that optimally stabilizes 
a targeted many-body system, which may be a crys- 
tal, disordered or quasicrystal structure, by incorporat- 
ing structural information that is not limited to the pair 



correlation function and generally accounts for complete 
configurational information. Unlike previous work, our 
primary interest is in the possible many-body structures 
that may be generated, some of which may include in- 
teresting but known structures, while others may repre- 
sent entirely new structural motifs. Moreover, the inverse 
methods described here can be employed to find targeted 
structures for metastable states as well as nonequilibrium 
configurations. 

In Section HH we define essential terms and briefly re- 
view basic concepts that are germane to the remainder 
of the article. Section IIIII describes inverse optimiza- 
tion methods for self-assembly of crystal ground states. 
In Section IIV1 we discuss recent applications of these 
methods that yield unusual crystal ground states with 
optimized, nondirectional interactions, including low- 
coordinated structures, chain-like arrays, and lattices of 
clusters. Section [V] describes and applies inverse opti- 
mization methods for self-assembly of disordered ground 
states. In Section I VII new duality relations for classi- 
cal ground states are reviewed and applied to some cases 
examined in the previous sections. Section [VIII discusses 
the pair-correlation- function realizability problem and in- 
verse optimization procedures to construct configurations 
with a given pair correlation function. In Section IVIII1 
we discuss inverse optimization methods to optimize in- 
teractions for targeted bulk properties and specific appli- 
cations. Finally, in Section HXl we suggest problems for 
future work and close with concluding remarks. 

II. BASIC DEFINITIONS AND CONCEPTS 

We consider a configuration of N identical interacting 
particles with coordinates = ri, r 2 , • • • , rjv in a region 
of volume V in (i-dimensional Euclidean space M. d . The 
coordinate of the ith particle generally embodies both 
its center-of-mass position and orientation as well as con- 
formation if required. In the absence of an external field, 
the classical iV-body potential <J>at can be decomposed 
into 2-body terms, 3-body etc., as follows: 

N N 

<f>N{r N ) = J2M riir ^ + ^ < Ps(? i ,T J) T h ) + ---. 

i<j i<j<k 

(2.1) 

Here ip n is the intrinsic n-body potential in excess to the 
contributions from <^2, ¥>3, • ■ • , <Pn—l- 

A given many-body structure is specified by the lo- 
cal density p(v), which can be expressed in terms of the 
particle coordinates as follows: 

N 

p(r) = J2^-r l ), (2.2) 

i=l 

where S(r) is the ^dimensional Dirac delta function. The 
N particle coordinates r N are statistically characterized 
by the ensemble (equilibrium or not) under considera- 
tion. The ensemble average of n products of the local 
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densities at n different positions yields the standard n- 
particle correlation functions^ For statistically homoge- 
neous systems in a volume V, these correlation functions 
are defined so that p n g n (r n ) is proportional to the prob- 
ability density for simultaneously finding n particles at 
locations r" = ri,ra, ...,r n within the system 7 -, where 
p = N/V is the number density. With this convention, 
each g n approaches unity when all particle positions be- 
come widely separated within V. Statistical homogeneity 
implies that g n is translationally invariant and therefore 
only depends on the relative displacements of the posi- 
tions with respect to some arbitrarily chosen origin of the 
system, i.e., 



fjn 



0n ( r i2, ri3, ■ ■ • ,r ln ), 



(2.3) 



where r ^ = r a — . 

The pair correlation function g% (r) is the one of pri- 
mary interest in this review. If the system is also ro- 
tationally invariant (statistically isotropic), then gi de- 
pends on the radial distance r = |r| only, i.e., 52(1") = 
92(f)- It is important to introduce the total correlation 
function h(r) = g^(r) — l, which, for a disordered system, 
decays to zero for large |r| sufficiently rapidly^ 

Such pair statistics can be inferred from radiation scat- 
tering experiments via the structure factor.— The struc- 
ture factor S(k), for an TV-particle system is related to 
the collective density variable 



N 



via the expression 



5(k) 



N 



(2.4) 



(2.5) 



where p(k) is the Fourier transform of p(r), defined by 
(|2.2|) , and i = Since the structure factor is propor- 

tional to the intensity of the scattered radiation, it is a 
nonnegative quantity for all k, i.e., 



S(k) > for all k. 



(2.6) 



This also mathematically follows from the nonnegative 
form (|2.5p . In the thermodynamic limit (N — > 00, V — > 
00 such that p is a fixed positive constant), the ensemble- 
averaged structure factor (omitting forward scattering) is 
defined by 



5(k) = l + ph(k), 



(2.7) 



where h(k) is the Fourier transform of the total correla- 
tion function h(r). 

The structure factor 5(k) provides a measure of the 
density fluctuations at a particular wave vector k. To 
see this important property quantitatively, consider the 
point pattern defined by the centers of particles in a 
many-body system at number density p. Let cr 2 (R) de- 
note the number variance of points contained within a 



d-dimensional spherical window of radius R in M. d . It can 
be shown 7 - 2 - that the number variance has the following 
real-space and Fourier-space representations: 



a 2 (R) = m (R) 



pvi{R) 



h(r)a(r; R) dr 



5(k)d(k; R)dk 



, (2.8) 



where v 1 (R) = n d / 2 R d /T(l + d/2) is the volume of the 
spherical window of radius R, a(r; R) is scaled intersec- 
tion volume, equal to the volume common to two spher- 
ical windows of radius R whose centers are separated by 
a distance r divided by V\(R), and a(k; R) is the corre- 
sponding Fourier transform of a(r; R). The scaled inter- 
section volume a(r; R) and its Fourier transform a(k; R) 
can be expressed explicitly in any dimension di 35 i 72 Thus, 
we see that the structure factor is directly related to the 
number variance at different wavelengths or, equivalently, 
for different window radii. In the limit of an infinitely 
large window, the relation above yields 



lim — y-t 

R^co pvi(R) 



= S(k = 0) = l + p h(r)dr. (2.9) 



Formula (|2.9[) applies whether the system is in equilib- 
rium or not. In the special case of an equilibrium sys- 
tem, it is well known that infinite-wavelength density 
fluctuations, as expressed by (|2.9|) . are proportional to 
the isothermal compressibility of the systemi 7 - 

For large R, it has been proved that a 2 (R) cannot grow 
more slowly than ^fR d ~ 1 or window surface area, where 
7 is a positive constant^ We note that point processes 
(translationally invariant or not) for which o~ 2 (R) grows 
more slowly than R d (i.e., window volume) for large R are 
examples of hyper uniform (or superhomogeneous) point 



patterns 



72,74 



Hyperuniformity implies that the structure 



factor S(k) has the following small k behavior: 



lim 5(k) 



0. 



(2.10) 



This classification includes all crystal structures^ point 
patterns associated with periodic and certain aperiodic 
tilings of space ) 72 ' 74 ! 75 ' 76 one-component plasmas ; 72 ' 74 
distribution of matter in the early Universe ) 77 ' 78 the 
ground state of superfluid helium, 79 maximally random 
jammed sphere packings^ and the ground states of spin- 
polarized fcrmions i 81 ' 82 



III. INVERSE METHODS FOR CRYSTAL 
GROUND STATES 

We recall that a classical ground-state configuration 
is one that minimizes the system potential en- 
ergy &n(t n . Our ability to identify ground states 
for a particular interaction is a highly challenging 
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problem. ; 83 ! 8 i 85 i 86 i 8 'i 88 not to mention the even more dif- 
ficult inverse problem of designing interactions to achieve 
targeted ground states. Here we describe recent progress 
on the latter problem. Because there is a vast (infinitely 
large) class of many-body potentials, we begin, for sim- 
plicity, by considering isotropic pairwise additive inter- 
actions, i.e., Eq. (1) reduces to the following form: 

N 

^ N (r N ) = ^(r ij ), (3.11) 

i<j 

where the pair potential ip(r) = f>2{r) is a radial func- 
tion, i.e., it depends on the radial distance r = |r|. Al- 
though realistic interactions that operate in soft matter 
systems can exhibit complicated many-particle charac- 
teristics, often a more economical description is sought 
that uses at most singlet and pair effective interactions 
that are density dependent to take advantage of the the- 
oretical and computational simplifications that result^ 
Therefore, our starting point of pairwise additivity (in 
the absence of an external field) is a practically useful 
approximation for colloids and polymers, for example. 

There are many open questions even for this sim- 
ple class of potentials. For instance, the limitations of 
isotropic pairwise additivity for producing target struc- 
tures are not fully known and can be probed using in- 
verse methods. We know that such interactions cannot 
produce thermodynamically stable chiral structures with 
a specified handedness; equal amounts of left-handed and 
right-handed structures would result. When is anisotropy 
in the potential required? An answer based on intuition 
from molecular systems would fail here. For example, the 
diamond crystal is thought to require directional interac- 
tions because such structures found in Nature result from 
covalent bonding. In fact, it has recently been shown 22 
that a diamond structure can be created from nondirec- 
tional interactions with strong short-range repulsions, as 
described in detail below. This structure has a special 
status in photonics research because a diamond crystal 
of dielectric spheres exhibits a photonic band gap across 
the Brillouin zone4^ 

Two inverse optimization schemes, called the "zero- 
temperature" and the "near-melting" schemes; 18 ' 19 have 
been devised for the purpose of designing interactions for 
targeted many-particle configurations. Specifically, the 
combination of these two optimization techniques lead 
to an A-body classical system with particles interacting 
via optimized potentials that has as its ground state (i.e., 
global energy minimum state) the corresponding target 
configuration in a specific volume (or density) range. Un- 
like previous attempts to solve this problem, this con- 
clusion is arrived at only after satisfying four important 
necessary criteria: 

1. lattice sums show that there is a positive pressure 
(or, equivalently density) range in which the given 
lattice is stable; 



2. all crystal normal mode frequencies are real at a 
specific density; 

3. defects (vacancies and interstitials) are shown to 
cost the system energy; 

4. and the system self-assembles in a molecular dy- 
namics or (Monte Carlo) simulation that starts 
above the freezing point and is slowly cooled. 

As concerns the last criterion, the system may start 
from an entirely random configuration or with a layer of 
fixed particles to promote epitaxial growth. Hence we 
make the important distinction here between homoge- 
neous and heterogeneous nucleation in self-assembly. It 
is of course a more stringent requirement that the de- 
sired lattice self-assemble from a random configuration 
(homogeneous nucleation). Note that sufficient criteria 
to ensure that a ground state is exactly achieved do not 
exist. 



A. Zero- Temperature Optimization Scheme 

In the zero-temperature optimization scheme, an op- 
timized pair potential for self-assembly of a partic- 
ular target configuration at a temperature of abso- 
lute zero is found by choosing a family of functions 
ip(r; Oq, a\, . . . , a n ), parametrized by the a^'s, and then 
finding the values of the parameters that lead to the most 
robust and defect-free self-assembly of the target crystal 
for a fixed density or, preferably, a range of densities (or, 
equivalently, a range of pressures). The objective func- 
tion is chosen so that the energetic stability of a given 
target crystal is maximized with respect to competitor 
lattices (chosen previously) subject to the condition that 
the target crystal is linearly mechanically stable. Me- 
chanical stability is ensured for a given potential by es- 
tablishing that that every phonon mode in the first Bril- 
louin zone is real. Thus, the structure is mechanically 
stable at zero temperature. However, this does not pre- 
clude other structures, periodic or otherwise, from being 
lower in energy than the targeted one. Therefore, the fi- 
nal outcome of the zero-temperature scheme becomes the 
initial potential function condition for the near-melting 
optimization procedure. 



B. Near- Melting Optimization Scheme 

From an initial parameterized potential (final outcome 
of the zero-temperature scheme) , the near- melting proce- 
dure optimizes the potential for self-assembly at a tem- 
perature near but below the crystal's melting point by 
suppressing nucleation of the liquid phase in molecular 
dynamics (MD) (or Monte Carlo) simulations. Specif- 
ically, simulations are repeatedly run at 80-95% of the 
melting temperature (the temperature is chosen such 
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that phase-transition fluctuations do not render the cal- 
culations inconsistent), each time calculating the Lindc- 
mann parameter, defined by 



Honeycomb Potential 



\ 



N ^ 



(3.12) 

where is the position of the i th particle (after an appro- 
priate amount of simulation time), is its initial posi- 
tion, and N is the number of particles. The parameter 
02 is then taken as the objective function for a simulated 
annealing calculation, and the ai are found such that O2 
is minimized. 

The ultimate criterion for self-assembly is the very 
strong condition that the targeted ground state be ob- 
served in a well-annealed molecular dynamics MD simu- 
lation starting from the liquid state. The system is slowly 
annealed to T = until the essentially defect-free target 
crystal results in reasonable computer time. Usually only 
a very few defects are found at the end of the simulation 
and its energy is checked to ensure that it is higher than 
that of the perfect crystal. 



IV. OPTIMIZED ISOTROPIC INTERACTIONS 
FOR LOW-COORDINATED CRYSTAL GROUND 
STATES 

Until recently, conventional wisdom presumed that 
low-coordinated crystal ground states require directional 
interactions. The aforementioned optimization schemes 
were tested initially to yield optimized isotropic (nondi- 
rectional) pair potentials that spontaneously yield the 
four-coordinated square lattice and three-coordinated 
honeycomb lattice as ground-state structures in two 
dimensions . 18 ! 19 The latter target choice is motivated by 
its three-dimensional analog, the diamond lattice. Figure 
[T] shows the optimized honeycomb potential and corre- 
sponding phonon spectra as well as annealed configura- 
tion at T = 0. It was found that as long as the salient 
features of the honeycomb potential are kept (two local 
minima at distance ratio v3j the first being positive and 
the second negative), self-assembly is unaffected by per- 
turbations in the potential; i.e., the potential is robust. 
This is an essential feature if this system is to be tested 
experimentally. We note that the functional form of the 
optimized "square-lattice" potential is simpler than that 
of the honeycomb crystal. 

In a separate work, the global phase diagram for the 
optimized "honeycomb" potential was determined. 913 The 
phase diagram was obtained from Hclmholtz free energies 
calculated using thermodynamic integration and Monte 
Carlo simulations. These results showed that the hon- 
eycomb crystal remains stable in the global phase dia- 
gram even after temperature effects are taken fully into 
account. Other stable phases in the phase diagram are 
high- and low-density triangular phases and a fluid phase. 




300 



250 




(c) 



FIG. 1: Honeycomb- crystal self-assembly as obtained in Ref. 
liH . (a) Optimized pair potential ip(r). Dimensionless energy 
and length units are defined by the axes for this potential, 
(b) Phonon spectrum (frequency squared) versus wave vector 
k for the optimized honeycomb-crystal potential at a specific 
area equal to 1.45. The acoustic and optical branches are 
shown, (c) Corresponding 500-particle annealed ground-state 
configuration. Although there are a few vacancies, these were 
"frozen in" during annealing due to finite time of the simu- 
lation. Such vacancies were shown to cost energy, indicating 
that the perfect honeycomb crystal is the true ground state. 
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FIG. 2: Unusual two-dimensional ground-state configura- 
tions generated from circularly symmetric pair potentials, as 
adapted from Ref. fTsL (a) Chain-like particle configurations, 
(b) Lattice of "simplex" clusters. 



bility than one would immediately think. It is also very 
possible that a much simpler isotropic potential could 
allow for a similar structure to assemble. 

These two-dimensional results were extended to the 
self-assembly of low-coordinated three-dimensional crys- 
tals with isotropic pair interactions, including the deter- 
mination of an optimized pair potential whose classical 
ground state is the simple cubic lattice and which is func- 
tionally simple enough to synthesize in the laboratory. 20 
The same investigation reported optimized isotropic po- 
tentials that yield the body-centered-cubic and simple 
hexagonal lattices (planes of triangular lattices stacked 
directly on top of one another), which provide other ex- 
amples of non-close-packed structures that can be assem- 
bled using only isotropic pair interactions. 

Optimized isotropic pair-interaction potentials with 
strongly repulsive cores have been obtained that cause 
the tetrahedrally coordinated diamond and wurtzite lat- 
tices to stabilize, as evidenced by lattice sums, phonon 
spectra, positive-energy defects, and self-assembly in 
classical molecular dynamics simulations. 22 Figure [3] de- 
picts one self-assembled diamond-crystal configuration 
shown from three different viewpoints. Finding such a 
potential via inverse methods is a highly nontrivial prob- 
lem, since the diamond crystal is extremely close in struc- 
ture to the tetrahedrally-coordinated wurtzite crystal in 
particular. Given the functional form of the potential, 
the pressure (or volume) was tuned very precisely to find 
a small stability range for the diamond structure, and 
under such conditions, simulations readily demonstrated 
its self-assembly. These results challenge conventional 
thinking that such open lattices can only be created via 
directional covalent interactions observed in nature and 
adds to our fundamental understanding of the nature of 
the solid state. 

Note that it has been shown that an isotropic pair po- 
tential that models star polymer systems has a region of 
phase stability that favors the diamond crystal. 10 How- 
ever, this potential, in contrast to the one reported in 
Ref. HH, possesses a soft core, which would be difficult to 
synthetically produce using colloids. 



No evidence of gas-liquid or liquid-liquid phase coexis- 
tence was found. 

In order to test the limitations of circularly symmetric 
potentials in two dimensions, pair potentials were devised 
that yielded structurally anisotropic chain-like arrays as 
well as lattices of compact clusters^ as shown in Fig. 
[2j These structures are reminiscent of "colloidal wires" 
and "colloidal clusters" found experimentally by the au- 
thors of Refs. [H and 0, respectively. Interestingly, we 
see that structural anisotropy (colloidal wires) can coun- 
tcrintuitively be achieved with isotropic interactions with 
the so-called "five-finger" potential^ This potential can- 
not be built in the lab with current technology, but it 
shows that isotropic potentials have perhaps more flexi- 



V. INVERSE OPTIMIZATION METHODS FOR 
DISORDERED GROUND STATES 

Collective density variables p(k) [cf. Q2.40 ] have proved 
to be useful tools in the study of static and dynamic phe- 
nomena occurring in many-body systemsi 91 ' 92 More re- 
cently, the collective-coordinate approach has been used 
to generate crystalline as well as noncrystalline classi- 
cal ground states for bounded or "soft" interactions us- 
ing numerical optimization techniques in two and three 
dimensions ! 17 i 27 ' 34 Soft interactions possess great impor- 
tance in soft-matter systems, such as colloids, microemul- 
sions, and polymers w*LL&2£idl 

We begin by briefly reviewing the basic description of 
the collective-coordinate approach for N interacting par- 
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(a)Vicw 1 




(b)View 2 



tides in a <i-dimensional cubic box of side length L and 
volume under periodic boundary conditions. The cor- 
responding infinite set of wave vectors is given by 



k=( 



2ixn\ 2irri2 



2im d 



(5.13) 



where the Hi (i — 1,2, ...,d) are positive or negative 
integers, or zero. The structure factor, defined by ()2.5|) . 
can be rewritten as follows: 



N 



1 



N 



C(k), 



where the real collective density variable C(k) is subject 
to the following constraints: 



~N(N-1) 



C(0) 
C(k) = C(-k) 

-^v<c(k) < 



(5.14) 
(5.15) 
(k^O). (5.16) 



The lower bound on C(k) arises because the structure 
factor 5(k) must be nonnegative. 

Let us assume that the total potential energy <5>at is 
pairwisc additive and therefore given by (|3.11[) . Suppose 
furthermore that the pair potential ip(r) has a Fourier 
transform ^(k): 

<^(k) = / tp(r) exp(ik • r)dr, (5-17) 
Jn 

tp(r) = fi" 1 ^(^(k)exp(-ik-r), (5.18) 




(c)View 3 

FIG. 3: Results of an MD simulation for 216 particles inter- 
acting via the optimized isotropic "diamond" potential show- 
ing self-assembly into a perfect diamond configuration, as in 
Ref. I22I . One configuration is shown from three different 
viewpoints, which clearly show that the result is the diamond 
crystal. 



where in the last expression the summation covers the 
entire set of k's. Then it is straightforward to show that 
the total potential energy for the iV-particle system can 
be exactly expressed in the following manner in terms of 
the real collective density variables 



<z> N (r N ) = n- 1 J2m)c(k). 



(5.19) 



Consider pair interactions whose transform y>(k) is non- 
negative radial function f(k) > (k = |k|) with compact 
support, i.e., 



where 



(p{k) = f(k)@(K - k), 



G(x) 




(5.20) 



(5.21) 



is the Heaviside step function. We see that if the C(k) 
is driven to its minimum value —N/2 for |k| < K, then 
that configuration must be a classical ground state of 
the system, the absolute minimum of <£>. Thus, density 
fluctuations for those k's such that |k| < K are com- 
pletely suppressed, i.e., the structure factor S(k) =0 for 
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|k| < K . Note that for the form (|5.20p . the corresponding 
real-space pair potential ip(r) will be an oscillating poten- 
tial. However, there are choices for f(k) one can make, 
especially for purposes of experimental realizability, that 
can appreciably dampen the amplitudes and range of the 
real-space interactions. 

Although the number of collective variables is infinite, 
the ./V-particle system possesses only dN configurational 
degrees of freedom, where d is the Euclidean space di- 
mension. Consequently, it is unreasonable to suppose 
(barring special circumstances) that generally all C(k)'s 
could be independently controlled. However, it is pos- 
sible, as illustrated below, to specify simultaneously a 
number of the collective variables equal to a significant 
fraction of dN. We denote by x the ratio of the con- 
strained degrees of freedom to the total number of de- 
grees of freedom. As \ increases to cover larger and larger 
numbers of the wave vectors, and consequently having an 
impact on a larger and larger fraction of the total degrees 
of freedom, the result for the classical ground state is far 
from obvious. It is clear that if x is a fraction of or- 
der unity, the ground state is periodic, which has been 
established^lM^Sa 

However, the more interesting cases involve disordered 
ground states (i.e., configurations that possess no long- 
range order), which arise for a range of x G [0, Xmax], pro- 
vided that Xmax is sufficiently smalh 17 ' 27 i 34 Our primary 
interest here are in the disordered, degenerate ground 
states that can be produced by the collective-density ap- 
proach. Such systems have the remarkable property of 
being able to self-assemble into one of the numerous de- 
generate disordered configurations when slowly cooled to 
absolute zero. 

For any given choice of N and K, the numerical pro- 
cedure utilizes a random number generator to create an 
initial configuration of the particles inside the hypercubic 
box. This starting point typically produces a large posi- 
tive value of the system potential energy <E> n . The next 
step involves use of an optimization procedure, such as 
the conjugate gradient method or a more sophisticated 
technique^ to seek a particle configuration that yields 
the absolute minimum value of $tv. 

This numerical optimization technique has been em- 
ployed to generate two-dimensional classical ground-state 
particle configurations with the simple transform choice 
f(k) = 1 in (|5.20p . i.e., the pure unit step function, which 
is zero for k > The resulting investigation distin- 

guished three structural regimes as the number of con- 
strained wave vectors is increased (i.e., as x is increased) 
- disordered, wavy crystalline, and crystalline regimes. 

The aforementioned collective-coordinate procedure 
has been generalized to those cases in which C(k) is con- 
strained to be some target value Cb(k) > for k e Q, 
where Q represents the finite set of k's for which a num- 
ber of the collective density variables can simultaneously 
be specified. Of course, each Co(k) must lie in the range 
specified by inequalities of (|5.16[) . Then consider the fol- 



lowing non-negative objective function: 

^(r N ) = E ^( k )P( k ) " C "(k)] 2 . (5.22) 

keQ 

If $jv is interpreted as a potential energy of interaction 
for the N point particles, then it can be shown that it 
represents intrinsic two-body, three-body and four-body 
interaction potentials operating in the system. If classical 
ground-state configurations for the N particles subject 
to that potential exist for which $n = 0, then those 
configurations necessarily attain the desired target values 
of the collective variables. 

This generalization of the collective coordinate ap- 
proach was applied in three dimensions^ In particular, 
multi-particle configurations were generated for which 
S(k) cx |k| Q , |k| < K, and a = 1, 2, 4, 6, 8, and 10. The 
case a = 1 is relevant for the Harrison-Zeldovich model of 
primordial density fluctuations of the early Universe ) 77 i 78 
superfluid helium^ maximally random jammed sphere 
packings^ and spin-polarized fermions i 81 ' 82 This analy- 
sis also provides specific examples of interaction poten- 
tials whose classical ground states for finite-sized systems 
are configurationally degenerate and disordered. 

Employing this collective-coordinate numerical opti- 
mization procedure, ground-state configurations of inter- 
acting particle systems in the first three space dimensions 
have been constructed so that the scattering of radiation 
exactly matches a prescribed pattern for a set of wave 
vectors i^l It is demonstrated that the constructed ground 
states are, countcrintuitively, disordered (i.e., possess no 
long-range order) in the infinite- volume limit. Three 
classes of configurations with unique radiation scatter- 
ing characteristics were studied: (i) "stealth" materials, 
which are transparent to incident radiation at certain 
wavelengths; (ii) "super-ideal" gases, which scatter radia- 
tion identically to that of an ensemble of ideal gas config- 
urations for a selected set of wave vectors; and (iii) "equi- 
luminous" materials, which scatter radiation equally in- 
tensely for a selected set of wave vectors. 

Although stealth materials and super-ideal gases are 
subsets of cqui-luminous materials, we use this term to 
refer to materials that scatter radiation more intensely 
relative to an ideal gas. These materials that scatter ra- 
diation much more intensely than an ideal gas for a set 
of wave vectors have enhanced density fluctuations and 
show local clustering similar to polymers and aggregating 
colloids i 94 ! 95 With the collective-coordinate inverse pro- 
cedure, the degree of clustering can be imposed by tuning 
the scattering characteristics for certain wavelengths. 

For purposes of illustration, disordered "stealth" con- 
figurations are depicted in two dimensions in Figure 2] 
for 168 particles for two selected values of x- At the low- 
est x considered, the configuration is seen not to have 
strong spatial correlations. At highest x value reported, 
the particles develop an exclusion shell about their cen- 
ters but the system still does not possess any long-range 
order. A system size study was carried out that revealed 
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FIG. 4: Stealth particle patterns of 168 particles in two di- 
mensions, as adapted from Ref. [U (a) \ ~ 0.04167, (b) 
X = 0.20238. Both systems are disordered but at higher \, 
particles tend to repel one another to a greater degree. The 
potential energy was minimized to within 10~ 17 of its global 
minimum. 



no long-range order when extrapolated to the infinite- 
volume limit. 



VI. DUALITY RELATIONS FOR CLASSICAL 
GROUND STATES 

The determination of the classical ground states of in- 
teracting many-particle systems (global minimum energy 
configurations) is an outstanding problem in condensed- 
matter physics and materials science . 83 ! 84 While classical 
ground states are readily produced by slowly freezing liq- 
uids in experiments and computer simulations, our the- 
oretical understanding of them is far from complete. 

Much of the progress to rigorously identify ground 
states for given interactions has been done for lattice 
models, primarily in one dimension^ 4 - The solutions in 
rf-dimensional Euclidean space M. d for d > 2 are consider- 
ably more challenging. For example, the ground state(s) 
for the well-known Lennard-Jones potential in R 2 or M 3 
are not known rigorously (although many computer sim- 
ulations support the conclusion that the hexagonal close- 
packed crystal is its ground state). 

Soft (bounded) interactions are easier to treat theo- 
retically as evidenced by recent progress in our under- 
standing of ground states of this class of potentials in R 2 
and R 3 i 17 i 27 i 34 i 86 i 87 Moreover, as we noted earlier, such 
interactions possess great importance in a variety of soft- 
matter systems 

Nonetheless, new theoretical tools are required to make 
further progress. Duality relations that link the energy of 
configurations associated with a class of soft pair poten- 
tials to the corresponding energy of the dual (Fourier- 
transformed) potential have recently been derived^ 
These duality relations enable one to use information 
about ground states of short-ranged potentials to draw 
new conclusions about the nature of the ground states 
of long-ranged potentials and vice versa. Among other 
results, they also have led to the identification of un- 
usual one-dimensional systems with ground-state "phase 
transitions" and can be employed to make computational 
searches for ground states more efficient. 

Before discussing these duality relations, which take 
the form of two theorems, we introduce some notation. 
Let U (r N ) be twice the total potential energy per particle 
in an iV-particle system with pairwise interactions, i.e., 



J7(r A ) 



hi 



(6.23) 



where tp(r) is a radial pair potential function and ry = 
|rj — A classical ground-state configuration is one 
that minimizes U(r N ). We consider those stable radial 
pair potentials tp(r) that are bounded and absolutely in- 
tegrable and call such functions admissible. Thus, the 
corresponding Fourier transform (p(k) in d dimensions at 
wave number k exists. We recall that in a Bravais lat- 
tice A, the space M. d can be geometrically divided into 
identical regions called fundamental cells, each of which 
contains one particle center^ We denote the reciprocal 
Bravais lattice by A. If the Bravais lattice A has density 
p, then its reciprocal lattice A has density p — p~ 1 (2ir)~ d . 
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Theorem 1. If an admissible pair potential ip(r) has a 
Bravais lattice A ground-state structure at number den- 
sity p, then we have the following duality relation for 
twice the minimized energy per particle U m i n : 

tp{r = 0) + = P0{k = 0) + p^W), (6.24) 

r£A keA 

where the prime on the sum denotes the zero vector 
should be omitted, A denotes the reciprocal Bravais lat- 
tice, and <p(k) is the dual pair potential, which auto- 
matically satisfies the stability condition, and therefore 
is admissible. Moreover, twice the minimized energy per 
particle U m in for any ground-state structure of the dual 
potential <p(k), is bounded from above by the correspond- 
ing real-space minimized quantity U m i n or, equivalcntly, 
the right side of If6.24p , i.e., 

Umin < U m in = P<f(k = 0) + p^ <P( fc )- (6.25) 

keA 

Whenever the reciprocal lattice A at reciprocal lattice 
density p = p~ 1 (2Tr)~ d is a ground state of (p{k), the 
inequality in !f6.25\) becomes an equality. On the other 
hand, if an admissible dual potential <p(k) has a Bravais 
lattice A at number density p, then 

U m in < U min = ptp(r = 0) + p~ 2 4 vK r )> (6.26) 

rGA 

where equality is achieved when the real-space ground 
state is the lattice A reciprocal to A. 

Whenever equality in relation (|6.25|) is achieved, then 
a ground state structure of the dual potential (p{k = r) 
evaluated at the real-space variable r is the Bravais lat- 
tice A at density p — p~ 1 (2n)~ d . Theorem 1 leads to 
another theorem (both of which are proved in Ref. l88l ) 
concerning phase coexistence. 

Theorem 2. Suppose that for admissible potentials 
there exists a range of densities over which the ground 
states are side by side coexistence of two distinct crys- 
tal structures whose parentage are two different Bravais 
lattices, then the strict inequalities in i6.25)) and i!6'.26'|) 
apply at any density in this density-coexistence interval. 

Note that the ground states referred to in Theorem 2 
are not only non-Bravais lattices, they are not even pe- 
riodic. The ground states are side-by-side coexistence of 
two crystal domains whose shapes and relative orienta- 
tions are complicated functions of p. 

On account of the "uncertainty principle" for Fourier 
pairs, a non-localized (long-ranged) potential ip(r) has a 
corresponding localized (compact) dual potential (p(k). 
Similarly, a localized (compact) potential <p(r) has a 
corresponding non-localized (long-ranged) dual potential 
(p(k). This property of Fourier pairs and the duality rela- 
tions of Theorem 1 enable one to use information about 



ground states of short-ranged potentials to draw new con- 
clusions about the nature of the ground states of long- 
ranged potentials and vice versa. In particular, three dif- 
ferent classes of admissible potential functions have been 
considered: (I) compactly supported functions (such as 
the ones employed in the collective-coordinate approach 
discussed in Section [Vj; (2) nonnegative functions; and 
(3) completely monotonic functions. 

For purposes of illustration, we discuss here in some 
detail, the application of Theorem 1 to the class po- 
tential functions that have been used in the collective- 
coordinate approach reviewed in Section[V] Recently, the 
ground states corresponding to a certain class of oscillat- 
ing real-space potentials <p(r) as defined by the family of 
Fourier transforms with compact support such that <p(k) 
is positive for < k < K and zero otherwise have been 
studiedJiSS Clearly, <p{k) is an admissible pair poten- 
tial. In Ref. HH, it was shown that in three dimensions 
the corresponding real-space potential ip(r), which oscil- 
lates about zero, has the body-centered cubic (bcc) lat- 
tice as its unique ground state at the real-space density 
p = l/(8-\/27r 3 ) (where we have taken K — 1). More- 
over, it was demonstrated^ that for densities greater 
than 1/(8v / 2tt 3 ), the ground states are degenerate such 
that the face-centered cubic (fee), simple hexagonal (sh), 
and simple cubic (sc) lattices are ground states at and 
above the respective densities l/(6\/37r 3 ), •\/3/(16v / 27r 3 ), 
and 1/(8V2 

Because all of the aforementioned ground states are 
Bravais lattices, the duality relation (I6.24|) can be ap- 
plied to infer the ground states of real-space potentials 
with compact support. Specifically, application of the 
duality theorem in R 3 and the results of Ref. |86| enables 
us to conclude that for the real-space potential ip(r) that 
is positive for < r < D and zero otherwise, the fee lat- 
tice (dual of the bcc lattice) is the unique ground state 
at the density y/2 and the ground states are degenerate 
such that the bcc, sh and sc lattices are ground states 
at and below the respective densities (3vo)/4, 2/v3, 
and l (taking D — l). Specific examples of such real- 
space potentials, for which the ground states are not 
rigorously known, include the "square-mound" potential 
[tp(r) = e > for < r < l and zero otherwise] and 
the "overlap" potential a(r;Z)/2)jS equal to the inter- 
section volume of two d-dimensional spheres of diameter 
D whose centers are separated by a distance r divided 
by the volume of a sphere (discussed in Section [TTJ) , and 
thus has support in the interval [0, D). Moreover, any 
structure, periodic or not, in which the nearest-neighbor 
distance is greater than unity is a ground state. 

Importantly, at densities corresponding to nearest- 
neighbor distances that are less than unity, the possi- 
ble ground-state structures is considerably more difficult 
to ascertain. For example, it has been argued in Ref. 
|43| (with good reason) that real-space potentials whose 
Fourier transforms oscillate about zero will exhibit poly- 
morphic crystal phases in which the particles that com- 
prise a cluster sit on top of each other. The square-mound 
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FIG. 5: (a) The localized square- mound potential [<p(r) = 
e — 1 for < r < 1 and zero otherwise] and overlap poten- 
tial [ip(r) = 1 — 3r/2 + r 3 /2 for < r < 1 and zero other- 
wise] in R 3 . (b) The delocalized dual square-mound potential 
(p{k) = tt 3/2 J 3/2 (fc)/(2fc) 3/2 multiplied by tt 3 /6 and dual over- 
lap potential (p{k) = 6n 2 J^ 2 (k/2)/k 3 . 



potential is a special case of this class of potentials and 
the fact that it is a simple piecewise constant function al- 
lows for a rigorous analysis of the clustered ground states 
for densities in which the nearest-neighbor distances are 
less than the distance at which the discontinuity in tp(r) 
occurs. 

The duality relations have also led to the identifica- 
tion of a one-dimensional system that exhibits an infi- 
nite number of "phase transitions" at T = from Bra- 
vais to non-Bravais lattices over the entire density range 
as well as a conjecture regarding the ground states of 
purely repulsive monotonic potentials^ Moreover, in- 
equalities (|6.25p and (|6.26[) provide a computational tool 
to estimate ground-state energies or eliminate candidate 
ground-state structures as obtained from annealing sim- 
ulations. 

The Gaussian potential is a special case of a purely re- 



pulsive monotonic potential, and is a useful interaction to 
model polymer systems ; 11 ! The phase diagram of such 
systems in various spatial dimensions has recently been 
investigated^ in order to understand the effect of dimen- 
sionality, apply the aforementioned duality relations, and 
to test a conjecture of Ref. [H concerning completely 
monotonic potentials. The Gaussian potential is an ex- 
ample of the class of potentials in which both the real- 
space and dual potentials are nonnegative functions. The 
authors of Ref. 4 - have argued that such systems display 
re-entrant melting with an upper freezing temperature. 

Elsewhere, corresponding duality relations for poten- 
tial functions that also include three-body and higher- 
order interactions will be derived^. 



VII. CONSTRUCTION OF CONFIGURATIONS 
WITH TARGET PAIR CORRELATIONS 

The subject of atomic and molecular distribution func- 
tions has enjoyed a long and rich history. However, not 
surprisingly for a scientific area so characterized by in- 
trinsic complexity, some deep problems of incomplete un- 
derstanding still persist. 

One such open question concerns realizability of a 
given candidate pair correlation function .92(1"), namely, 
whether it actually represents the pair correlation of some 
many-particle configuration at number density p > 0. 
This is called the realizability problem i 30 ' 72 Several nec- 
essary conditions that must be satisfied by the candidate 
are known, including nonnegativity of g% (r) and its asso- 
ciated structure factor S(k), as well as constraints on im- 
plied local density fluctuations^ It has recently come to 
light that a positive at a positive p must satisfy an un- 
countable number of necessary and sufficient conditions 
for it to correspond to a realizable point process i 100 ' 101 
However, these conditions are very difficult (or, more 
likely, impossible) to check for arbitrary dimension. In 
other words, given p and <?2, it is difficult to ascertain if 
there are some higher-order functions 33 , 54 , . . . for which 
these one- and two-particle correlation functions hold. 

To shed light on the realizability problem, a simple one- 
dimensional lattice model, with single-site occupancy, 
and nearest-neighbor exclusion has been investigated. 32 
The following results were obtained: (a) pair correla- 
tion realizability over a nonzero density range, (b) vi- 
olation of the Kirkwood superposition approximation for 
(73, and (c) inappropriateness of the so-called "reverse 
Monte Carlo" method that uses a candidate pair corre- 
lation function as a means to suggest typical many-body 
configurations. Note that Chayes and Chayes 1 ^ 2 - proved 
that for any pair correlation function (meeting mild con- 
ditions) that is derivable from an iV-body Hamiltonian, 
there always exists a unique "effective" two-body po- 
tential that produces the same pair correlation function 
(but generally not the higher-body correlation function 
S3 7 <?4 j <?5 1 etc.). This theorem has been successfully ap- 
plied to polymer solutions to obtain effective pair inter- 
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actions from g ^ 103 ' 104 - 

Elsewhere, so-called iso-g2 processes were studied in 
the equilibrium regime. These consist of a sequence 
of equilibrium many-body systems that have different 
number densities but share, at a given temperature, the 
same "target" pair correlation function. In other words, 
in these processes, density-dependent interactions iden- 
tically cancel the usual density variation of many-body 
pair correlation functions ! 28 ' 29 ! 33 Target pair correlation 
functions studied include the unit step function as well 
as the zero-density limit of the square- well potential (for 
which .92(f) = exp[— /3<p(r)]). Formal density expansions 
for effective pair potentials were derived with this iso- 
g2 property showing how successive terms in that ex- 
pansion can be determined iteratively. Explicit results 
through second density order have been obtained for two 
types of "target" pair correlation functions, and the con- 
ditions under which realizability can be attained were 
explored.— 

In order to explore and gain insight into the basic sta- 
tistical geometric features of random sphere packings, the 
notion of a (72-invariant process was introduced. 30 A #2- 
invariant process is one in which a given nonnegative pair 
correlation (72 (r) function remains invariant as density 
varies for all r over the range of densities 



< p < p* 



(7.27) 



The terminal density is the maximum achievable den- 
sity for the t^-invariant process subject to satisfaction of 
the known necessary conditions on the pair correlation 
function. The determination of the terminal density for 
various forms of g-i that putatively correspond to a sphere 
packing has been solved using numerical and analytical 
optimization techniques i 30 ' 35 ' 36 ' 37 

To test whether such gi 's at terminal density are in- 
deed realizable by sphere packings, stochastic optimiza- 
tion techniques originally, developed to construct ma- 
terial microstructures with targeted lower-order corre- 
lation functions ; 105 ' 106 ' 107 were employed. 31,34 In a con- 
struction algorithm, an initial configuration of particles 
evolves such that the final configuration possesses a set of 
targeted correlation functions up to some "cut-off" dis- 
tances. This is done by choosing the objective function to 
be a "squared error" involving the set of targeted correla- 
tion functions. The evolving configurations are induced 
by minimizing this objective function via a stochastic op- 
timization procedure. 

For example, for the case of ci-dimensional packing of 
congruent spheres of diameter D in which the pair cor- 
relation function is taken to be 



52 M = 



z 



psi(r) 



6(r-D)+Q(r-D), 



(7.28) 



where Z represents the average contact value per sphere 
and si(r) = dn^ 2 ^' 1 + d/2) is the surface area of 
a d-dimensional sphere, 3 & it was found that the termi- 
nal packing fraction (fraction of space covered by the 
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FIG. 6: (a) Graph of the target pair correlation function (72 (r): 
Dirac 8 function plus a step function, (b) A two-dimensional 
configuration of 500 particles that realizes this targeted form 
for g2(r) up to a dimensionless distance of r/D = 2.5, as 
adapted from Ref. 03. The configuration consists of only 
dimers at the terminal packing fraction 0* = 0.5 with an 
average contact value Z = 1.0. 



spheres) and the associated average contact number 
are given by 



d + 2 



Z M 



d 
2' 



(7.29) 



Numerical evidence suggests that such a pair correlation 
is achieved by a single sphere packing configuration for 
any d > 2.— Such a pair correlation function Figure [6] 
shows a realization of such a packing in two dimensions. 
Of course, in any simulation, pair distances must binned 
and sampled up to some cut-off distance. Note that for 
a sufficiently large system, the targeted correlation for 
a single configuration approaches that of one obtained 
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from an ensemble of configurations by ergodicity. 

Because the realizability problem is far from being 
solved, it remains an active area of research. For exam- 
ple, it has been conjectured that any radial, nonnegative 
pair correlation function characterized by a hard-core, 
which decays sufficiently rapidly to unity, is realizable by 
a translationally invariant disordered sphere packing in 
d-dimensional Euclidean space for asymptotically large d 
if and only if S(k) > 0. 108 Although there is mounting 
evidence to support this conjecture ; 35 ' 37 ' 108 a proof of it 
is a great challenge. 



VIII. DESIGNING ISOTROPIC PAIR 
POTENTIALS FOR TARGETED BULK 
PROPERTIES 



teractions. (Note that NTE has been shown to occur in a 
two-dimensional fluid with isotropic interactions ! 112 ' 113 ) 
Moreover, it was established that a sufficient condition 
for a potential to give rise to a system with NTE behavior 
is that it exhibits a softened interior core within a basin 
of attraction (as depicted schematically in part (a) of Fig. 
|7J|. Using an optimization procedure to find a potential 
that yields a strong NTE effect and constant-pressure 
Monte Carlo simulations, it was shown that as the tem- 
perature was increased, the "softened interior core" po- 
tential [part (b) of Fig. [7], the system exhibited negative, 
zero, and then positive thermal expansion before melting 
(in both two and three dimensions). 



Inverse methods have been recently devised to optimize 
interactions of many-particle systems to achieve targeted 
novel bulk properties. To illustrate the interesting possi- 
bilities, we discuss three specific target examples in some 
detail: the thermal expansion coefficients and Poisson's 
ratio. 



A. Thermal Expansion Coefficients 

Control of thermal expansion properties of materials 
is of technological importance due the need for struc- 
tures to withstand ambient temperature variations. In 
the technological realm, materials with zero thermal ex- 
pansion (those that do not expand or contract upon heat- 
ing) can aid in the longevity of space structures, bridges 
and piping systems i 50 ' 109 Materials with very large ther- 
mal expansion coefficients could function as actuators, 
and those with negative thermal expansion coefficients 
may be of use as thermal fasteners^ 8 - 

Negative thermal expansion (NTE) behavior, a well- 
known but unusual phenomenon in many-particle sys- 
tems, has been observed only in multi-component mate- 
rials with open unit cell structures in which the bonding 
of component particles is highly directional. Perhaps the 
most common example of a solid exhibiting NTE is that 
of ice, which contracts upon melting into liquid water J^- 
Another example of a material that undergoes NTE is 
zirconium tungstate, ZrWiO%, which exhibits this be- 
havior for an extremely large temperature range, namely 
0.3K through 1050^^ 

An isotropic interaction potential has been optimized 
that gives rise to negative thermal expansion (NTE) 
behavior in equilibrium many-particle systems in the 
solid state in both two and three dimensions over a 
wide temperature and pressure range (including zero 
pressure) i^ 3 - Although such anomalous behavior is well- 
known in materials with directional interactions (e.g., 
zirconium tungstate), this is the first time that NTE be- 
havior has been established to occur in the solid state of 
single-component many-particle systems for isotropic in- 
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FIG. 7: (a) Schematic depiction of an isotropic pair potential 
(scaled by the well depth e) with a softened interior in its 
basin of attraction following Ref. [U Thermal fluctuations 
cause the average nearest-neighbor distance to decrease, re- 
sulting in an overall contraction of the system upon heating, 
(b) The optimized "softened interior core" (SIC) potential, as 
adapted from Ref. [23l . has NTE behavior over a wide range 
of temperatures. 
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B. Poisson's Ratio 

Another interesting target bulk property is the Pois 
son's ratio v. In particular, it is desired to optimize in 
teractions to achieve negative Poisson's ratio (NPR), th 
so-called "auxetic" material a 114 ' 115 . When such materi 
als are stretched in a particular direction, they expand h 
an orthogonal direction. Auxetic behavior is a counter 
intuitive material property that has been observed onl; 
in a handful of elastically isotropic materials that oftei 
have intricate structures and characteristic lengths mud 
larger than an atomic bond length, such as foams 5 - 1 - an< 
other cellular materials i 52 i 53 ' 54 i 56 Auxetic materials hav 
a great deal of technological potential; for example, the; 
can be used as strain amplifiers. 115 If auxetic materi 
als are used as a matrix in the manufacture of minia- 
ture sensors based on piezoceramic composites, the range 
of operating frequencies of a piezoelectric transducer is 
widened and the sensitivity of the device is increased; 1 ^ 
They can also be used as mechanical components of mi- 
croelectromechanical systems, and as transducing struc- 
tures, shock absorbers and fasteners^ 

It has been recently found that under tension (i.e., neg- 
ative pressure), many-body two- and three-dimensional 
systems with isotropic two-body interaction potentials 
can have a negative Poisson's ratio in the crystal phase 
as long as certain linear equalities and inequalities in- 
volving the interaction potential ip(r) are satisfied^ This 
is an unexpected result, since it describes an inherently 
anisotropic behavior that arises from isotropic interac- 
tions; indeed, most previously discovered auxetic mate- 
rials exhibit complex, carefully designed anisotropic in- 
teractions. This can be shown to be the case at zero 
temperature for the elastically isotropic triangular lat- 
tice in two dimensions, and for the fee lattice in three 
dimensions, which, surprisingly, can also be made to be 
elastically isotropic. One can show that in the former 
case, the simple Lennard- Jones potential can give rise to 
auxetic behavior (see Fig. [8]). In the three-dimensional 
case, auxetic behavior is exhibited even when the elastic 
constants are constrained such that the material is elas- 
tically isotropic. Finding auxetic behavior over a wide 
range in temperature and pressure is a challenging opti- 
mization problem that has yet to be addressed. 

This analysis suggests that auxetic behavior only oc- 
curs in crystals under the nonequilibrium condition of 
negative pressures when the system contains only pair 
interactions and is elastically isotropic. Such auxetic 
materials may potentially be experimentally produced 
using synthetic techniques that rely on kinetic effects; 
examples include tempered glasSfiii and even colloidal 
crystalsi 1 ^ However, a three-body potential has been de- 
vised that yields NPR behavior in close-packed two- and 
three- dimensional lattices by construction at zero tem- 
perature and positive pressured In order to produce this 
behavior, the potential has a built-in energy cost associ- 
ated with deforming the equilateral triangles in the two- 
dimensional triangular lattice and the three-dimensional 
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FIG. 8: Region of lattice constants (indicated by the rect- 
angular box) for which the Poisson's ratio is negative in a 
triangular lattice, using the Lennard- Jones interaction poten- 
tial, (pLj, as adapted from Ref. [26| Pressure is positive to the 
left of the dotted line and negative to the right; thus, auxetic 
behavior only occurs at negative pressure. To the right of the 
rectangular box, the lattice becomes unstable. 

close-packed lattices. The interested reader is referred to 
Ref. [25| for the explicit form of this three-body potential. 

IX. FUTURE WORK AND CONCLUSIONS 

In this section, we discuss future directions and close 
with concluding remarks. 

A. Interaction Potentials for Targeted 
Configurations at Positive Temperature 

Most of the inverse techniques reviewed here were di- 
rected toward obtaining ground state (T = 0) struc- 
tures. However, the same methods can be extended 
to treat many-particle configurations at positive tem- 
perature. For example, an ability to control the for- 
mation of point, line, and planar defects of crystals 
under various growth conditions at positive tempera- 
ture is highly desirable. The required interactions to 
achieve representative amorphous target structures, in- 
cluding equilibrium liquids at positive temperature and 
low-temperature glasses, is another interesting applica- 
tion. 



B. Interaction Potentials for Targeted 
Multicomponent Systems 

It is straightforward to extend the zero-temperature 
and near-melting optimization scheme s 18 ' 19 to multicom- 
ponent systems. The parameter space, which now in- 
cludes species composition and effective particle size ra- 
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tios, becomes much larger than the single-component in- 
stance, and therefore one must be careful in selecting 
the family of potential functions that must be optimized 
as well as the target structures. In order to make the 
search manageable, one could limit the choice of poten- 
tial functions to those that are consistent with interpar- 
ticle interactions found in colloidal systems. In addition 
to hard-sphere-like interactions, these include long-range 
repulsive, short-range attractive and averaged dipolar in- 
teractions. It has recently been shown that the electro- 
static interaction between oppositely charged particles, 
which arc long-range attractive interactions, can result 
in a rich class of stable ionic colloidal crystals^ as illus- 
trated Figure El Motivated by this remarkable investi- 
gation, one can imagine optimizing a family of potential 
functions based on such interactions that target an even 
broader class of crystal structures. 




FIG. 9: Fheoretically predicted stable binary crystals of op- 
positely charged colloids with different stoichiometries, as ob- 
tained from Ref. [l4| (with the permission of the authors). 



C. Anisotropic Interactions 

We have seen that there exists nontrivial families of ra- 
dial pair potentials for which interesting targeted struc- 
tures are the stable low-temperature forms. Conse- 
quently, it was not necessary in principle to call upon 
angle-dependent or non-additive interactions to form 
such nonconventional lattices. However, anisotropic pair 
interactions offer greater flexibility to achieve targeted 
structures and therefore provides a new direction to ap- 
ply our inverse methods. 

Recently, a new generation of colloidal particles with 
chemically or physically patterned surfaces has been de- 
signed and synthesized in the attempt to manipulate the 
valency of the colloidal particles.— ' 119 ' 120 This synthesis 
effort aims to generate "superatoms" (i.e., atoms at the 
nano and microscopic length scales) in order to repro- 
duce and extend traditional collective molecular behav- 



ior to larger length scales; thus, opening the new eld of 
"supra-particle" colloidal physics. 

One simple way to model such interactions is via 
"patchy" particles, i.e., particles with discrete, attractive 
interaction sites at prescribed locations on the particle 
surface. Molecular simulations have been carried out to 
investigate the self-assembly of patchy particles i 12 ' 16 ' 121 
Chains, sheets, rings, icosahedra, square pyramids, tetra- 
hedra, and twisted and staircase structures have been 
obtained through suitable design of the surface pattern 
of patches. Patchy particles represent a new class of 
building blocks for the fabrication of colloids with unique 
structural characteristics. 

Thus, it would be highly desirable to optimize patchy 
particles interactions to achieve low-coordinated crys- 
tal structures, amorphous structures, and quasicrystals. 
Again, this can be accomplished by appropriate simple 
extensions of the zero-temperature and near-melting op- 
timization schemes that were originally implemented for 
isotropic interactions. 



D. Inverse Optimization Methods for Novel 
Targeted Bulk Properties 

A full-blown and general optimization scheme that can 
be used to find optimized interactions over a large fam- 
ily of potential functions for a given set of bulk prop- 
erties over a wide range of conditions has yet to be de- 
vised. For example, optimizing interactions in a many- 
particle system so that it exhibits auxetic behavior over 
a wide range in temperature and pressure is a challeng- 
ing problem. One path toward the general goal is to 
formulate a methodology that incorporates a set of bulk 
properties in the objective function in the same spirit as 
has been done for topology optimization of composite 
materials) 48 ' 59 ' 61 ' 109 but in a molecular dynamics sim- 
ulation. Specifically, the objective function can either 
be the bulk property itself (which is extremized) or a 
squared "error" function involving a targeted bulk prop- 
erty (which is minimized) during a molecular dynamics 
simulation. The simulation would start from some initial 
configuration and randomly distributed velocities for a 
initial parameterized potential. At fixed time intervals, 
the objective function would be computed and then the 
parameters of the potential updated according to some 
optimization routine (e.g., simulated annealing). This 
procedure would then be iterated until the objective func- 
tion is extremized. 

An intriguing set of target materials are those that 
exhibit "inverse melting."— Inverse melting is a first- 
order phase transition involving the crystal and liquid, 
but with a reversal from conventional melting in that ad- 
dition of heat to the liquid, at constant pressure, causes 
that liquid to freeze into a crystalline solid. As a re- 
sult of this reversal, the crystal has higher entropy than 
the isotropic liquid with which it coexists. This is a 
rare phenomenon, but real- world examples exist. For 
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example, the transition itself forms the basis of the zone- 
refining method for purification. 122 Inverse melting has 
been studied as a "forward" problem using the Gaussian 
core potential model. 67 However, devising optimized in- 
teractions to make such unusual macroscopic behavior 
as robust as possible over a wide range of conditions has 
heretofore not been considered. 



E. Toward Experimentally Realistic Interactions 

An important component of future research should be 
the development of robust potentials (even if not op- 
timal) for targeted structures and bulk properties that 
can be synthesized experimentally with colloids or other 
soft-matter systems. It is clear that there is wide class 
of target structures and bulk properties that can be 
achieved with pairwise additive potentials, both isotropic 
and anisotropic. In future research, it will be highly de- 
sirable to determine, when possible, the pair interactions 
that can be either be synthesized experimentally with 
colloids using current technology (e.g., depletion, screen- 
ing length, dipolar interactions, etc.) or can be done so 
in the near future. The latter could serve as a challenge 
to experimentalists. 

Real interactions in many-particle materials at nondi- 
lute concentrations are necessarily nonadditive, i.e., 
intrinsic three-body and higher-order interactions be- 
yond pair interactions [explicitly given in Eq. (1)] are 
inevitable * 89 ' 104 Thus, it is crucial to determine how the 
effective pair potentials that result from the inverse ap- 
proach correspond to the many-body interactions that 
arise in actual colloidal systems. This important problem 
has received little attention in the literature. It has been 
shown that effective pair interactions that approximate 
nonadditive potentials are in fact density dependent and 
hence one must be careful in carrying out the resulting 
statistical mechanics*^, Guided by experiments, one can 
determine the real two-body and three-body interactions 
that together mimic the effective pair potential required 
to achieve the targeted many-particle configurations us- 
ing both theoretical techniques and molecular dynamics 
simulations. This will require continual feedback between 
theory and experiment. 



F. Incorporating Dynamics 

The dominant theme of this review article concerned 
the determination of potentials that spontaneously create 
target structures under equilibrium or near-equilibrium 
circumstances. A conjugate kinetic problem also exists, 
in which selection among alternative irreversible scenar- 
ios (involving distinct dynamical evolutions) itself be- 
comes a tool for selecting among alternative structural 
outcomes. The full potential of self-assembly to control 
and manipulate the structure of materials at the micro- 
scopic and nanoscopic level cannot be realized without 



a deeper understanding of nonequilibrium processes at 
those length scales. For example, a recently developed 
model demonstrates this point by showing how the ir- 
reversible collisions in particle suspensions that gener- 
ally produce diffusive chaotic dynamics can also cause 
the system to self-organize to avoid future collisions 
This can lead to a self-organized non-fluctuating quies- 
cent state, with a dynamical phase transition separating 
it from fluctuating diffusing states. This investigation 
and many other nonequilibrium studies, too numerous 
to list here, provide exciting glimpses into the future of 
self-assembly Inverse optimization techniques that ex- 
ploit the dynamics of many-particle systems to achieve 
self-assembly has yet to be developed and should offer 
greater flexibility for novel material design. 

G. Conclusions 

Although in their infancy, the inverse approaches re- 
viewed here have already shown a capability for con- 
trolling self-assembly to an exquisite degree. Indeed, 
future applications could revolutionize the manner in 
which materials are designed and fabricated, especially 
if there is continual feedback between theory and ex- 
periment. There are recent examples in which output 
from material optimization studies have been combined 
with experiments to produce novel materials or material 
components . 124 ' 125 ' 126 ' 127 These inverse methods have led 
to a deeper fundamental understanding of the mathemat- 
ical relationship between the collective structural behav- 
ior of many-body systems and their interactions. For ex- 
ample, we have seen that low-coordinated crystal struc- 
tures, chain-like arrays, and layered structures do not re- 
quire directional interactions for self-assembl y 18 ' 19 ' 20 ' 22 
Although soft matter with some of the interactions re- 
viewed in this article cannot be synthesized with cur- 
rent technology, other optimized interactions described 
here that yield either novel structures or bulk proper- 
ties are rather standard or could easily made in the 
laborator y 20 ' 23 ' 26 . For practical purposes, it will be im- 
portant that future research be directed toward produc- 
ing optimized interactions with the constraint that they 
are experimentally achievable. We envision being able 
to "tailor" potentials that result in novel materials with 
varying degrees of disorder, thus extending the tradi- 
tional idea of self-assembly to incorporate not only crys- 
tals but amorphous and quasicrystal structures. 
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